Polarization forces in water deduced from single molecule data 
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Intermolecular polarization interactions in water are determined using a minimal atomic multipole 
model constructed with distributed polarizabilities. Hydrogen bonding and other properties of 
water- water interactions are reproduced to fine detail by only three multipoles /j,h, mo, and 8o and 
two polarizabilities cto and an, which characterize a single water molecule and are deduced from 
single molecule data. 
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Understanding polarization forces is crucial in many 
molecular systems such as molecular clusters, liquids, 
or solids, specifically those containing polar and polar- 
izable molecules. Polarization effects in water are partic- 
ularly strong, as can be judged by the enhancement of the 
molecular dipole from 1.855 D for an isolated molecule to 
2.6 — 3.2 D in condensed state. 1 ' 2 Water is a very basic 
substance. 3 It is a fascinating object to study because of 
its singular properties, its significance in biological sys- 
tems, and because it is a classic example of hydrogen 
bonding. 4 Hydrogen bonding, which itself is one of the 
key elements of the functioning of life, is largely a po- 
larization effect. Unfortunately, no commonly accepted 
model describes it simply and accurately at the same 
time. Here we show that application of recent rules for 
minimal atomic multipoles 5 combined with the notion of 
distributed polarizabilities lead straightforwardly, with- 
out further intervention, to a very transparent model for 
polarization forces in water. Hydrogen bonding and other 
properties of water- water interactions are reproduced to 
fine detail with only three atomic multipoles and two po- 
larizabilities, whose values are deduced based on single 
molecule data. 

Intermolecular potential for water has been extensively 
studied, with about 150 models introduced since 1930s, 
indicating difficulties in this area. 6 Recent accurate pa- 
rameterizations involving several tens of parameters are 
available based on tuning to rich vibration-rotation- 
tunneling (VRT) spectra, 7,8 or to high-level quantum- 
chemical calculations, 9 or both. 10 Some models are based 
on molecular multipole moments and require high-order 
multipoles. 11 Following seminal work by Rahman and 
Stillinger, 12 many empirical models involve distributed 
charges. 13-22 Most of the force fields use static charges 
thus ignoring or averaging the polarization effects, while 
other models incorporate polarizabilities explicitly. 16-22 
Work 17 first introduced molecular polarizability of water 
distributed over atomic sites. 

It has been recently recognized that hydrogens need 
not be assigned charges in distributed charge models. 5 
The hydrogen's sole electron participates in the chemical 
bond and is not centered at the proton. Therefore, hy- 
drogen is best described by an atomic dipole placed at 
the proton and directed along the bond. Assigning both 



charge and dipole causes redundancy and leads to un- 
physical results. This rule is an integral part of the mini- 
mal atomic multipole expansion (MAME), 5 which elimi- 
nates the redundancies by a careful choice of the minimal 
set of atomic multipoles based on the Lewis structure of 
the molecule. 

MAME rules lead to the following expression for the 
electrostatic potential of a single water molecule: 
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Since protons have no charge, neutrality allows no charge 
on the oxygen either. The dipole fio and quadrupole 
9o describe the two lone pairs on oxygen. 5 Origin is at 
the oxygen, ri j2 are the positions of protons, n, 2 = I, 
n = (ri+r 2 )/|ri+r 2 | is the unit vector along the symme- 
try axis, and ni. 2 are unit vectors in the directions of lone 
pairs (Fig. 1). Experimental geometry has I = 0.9572 A 
and a nearly tetrahedral bond angle (3 = 104.52° be- 
tween ri and r 2 . 23 We take and n 2 to be at the tetra- 
hedral angle f3 = 109. 47°. 19 Significant deviation from 
this value leads to a dramatic deterioration of accuracy 
of Eq. (1) as seen in the inset. 




p' = 109.47 



Fig. 1 Geometry of a single water molecule. Vectors ni 
and ri2 point in the direction of the lone pairs on oxygen. 

The inset shows an effect of varying /3 on the accuracy 5 of 

(1). /in, fio and 9o are re-optimised for every j3 . The vertical 
bar marks the perfect tetrahedral angle. 
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Our goal is to extend the static model to describe the 
field induced by a polarized molecule subject to external 
fields. In doing so we again keep only the minimal set 
of multipoles to avoid redundancies. Charge redistribu- 
tion part of the molecular polarizability 24 vanishes for 
water due to the absence of charged sites. Thus we as- 
sign polarizabilities to individual nuclei in such a way as 
to reproduce experimental molecular polarizability. The 
smallest component is a yy = 1.4146(3) A 3 normal to the 
molecular plane, the next is a zz = 1.4679(13) A 3 along 
the dipole moment, and the largest is a xx = 1.5284(13) 
A 3 in the longest dimension. 25 

Atomic polarizabilities reflect the local atomic environ- 
ments and need not necessarily be isotropic. Tetrahedral 
coordination of oxygen suggests to assign it an isotropic 
polarizability ao- For hydrogens the polarizability o>h 
along the OH bond may differ from the polarizability a± 
normal to it. To deduce ao, an, and a± we express the 
molecular polarizability, 

a xx = a + 2a H sin 2 0/2 + 2a± cos 2 0/2, 
oiyy = a + 2a± , 

a zz = a + 2a H cos 2 (3/2 + 2a ± sin 2 0/2. (2) 

In a surprise twist, the determinant of this linear system 
is identically zero. Equations (2) are therefore dependent 
and possess a solution only if the quantity 

a xx cos 2 0/2 + a yy (2 sin 2 0/2 - 1) - a zz sin 2 0/2 (3) 

is zero. Thus, the model is adequate if the relation 
holds between the molecular polarizability components. 
Substituting the experimental values into (3) we get 
0.0093A 3 , which is indeed close to zero. Two indepen- 
dent equations suggest that one of the atomic polariz- 
abilities can be safely omitted. The natural choice is to 
set a±_ — 0, implying that the dipole moments on protons 
can change their value, but not direction. Solving (2) we 
get 

a Q = a yy — 1.4146 A 3 , and 

% = (a xx + a zz )/2 - u yy = 0.0836 A 3 . (4) 

Thus, the bulk of molecular polarizability comes from the 
oxygen, which is consistent with its atomic size, while 
the small polarizabilities on the protons account for the 
(small) anisotropy of the molecular polarizability tensor. 

Three gas-phase multipoles from a density functional 
calculation, fiH = 0.675 D, no — 1-033 D, and &o — 
1.260 DA 5 result in the molecular dipole fi = 1.854 D 
and the quadrupole components O = <d xx — Q yy = 4.973 
DA, Z2 = 0.142 DA. 36 These should be compared to 
experimental data, 26 ' 27 n = 1.8546(6) D, 6 = 5.126(25) 
DA, and e zz = 0.113(27) DA. 

We again adjust the three atomic multipoles to satisfy 
the three experimental values precisely to avoid any com- 
putational input. The molecular dipole and quadrupole 
are expressed in terms of the atomic multipoles as 



H = no + 2fi H cos/3/2, 
6 = 6lfi H sin 2 0/2 + W sin 2 0/2, 
Qzz = 2^ H (3cos 2 /5/2 - 1) -6» (3 cos 2 /?'/2- 1) (5) 

In practice, we face here an almost identical problem, in 
that the determinant of (5) is small. It becomes zero 
when an ideal tetrahedral angle is substituted for 0. A 
relation similar to (3) in this case reads simply <3 ZZ = 0. 
Actual Q zz is indeed small, but not zero, and deviates 
noticeably from 109.47°. Nevertheless, smallness of the 
determinant indicates that the finite accuracy data can 
be satisfied by a range of atomic multipoles, and so the 
third equation in (5) cannot be used reliably. 

Thus, we use the first two equations to express no and 
80 in terms of nn , which guarantees to reproduce experi- 
mental n and O, while keeping reasonable Q zz . The DFT 
value hh = 0.675 D yields no = 1-029 D and 6 Q = 1.352 
DA, with Q zz = 0.160 DA. The model is thus completely 
defined and readily yields the polarization energy Ep for 
the water dimer, trimer and larger clusters. 24 

Water clusters from dimers on up have been ex- 
tensively studied with both experiment 7,28-30 and 
theory. 7,8,31-33 Six-dimensional adiabatic energy surface 
of the dimer has 8 equivalent minima 34 split in a com- 
plex fashion by zero-point tunneling motion. Softness of 
the pair potential requires care when relating it to the 
experimental observables. 7 





i 1 1 

\ & -<z 

j E P (e d ) \ 


... \ s~ 








R 1 






/ \ E p (9 a ) \ . 






" 6 


HI ... _ 

■ III 


1 ■ 


■ 



-120 -90 -60 -30 30 60 90 

Fig. 2 Polarization energy (Coulomb + induction), 
Ep(6 a , 0d) for water dimer. Roo = 2.977A, and one of the an- 
gles is fixed at the minimum value, as the other one is varied. 
Ep is calculated for the defined system of atomic multipoles 
and polarizabilities in the standard manner, 24 by computing 
fields of all multipoles of one molecule exerted on the multi- 
poles of another molecule, and solving for self-consistency. 

Equilibrium hydrogen bonded configuration has a sym- 
metry plane (Fig. 2, inset) and is characterized by the 
oxygen-oxygen distance Roo , the donor angle 9d and the 
acceptor angle 9 a . The hydrogen bond forms when the 
donor proton points against one of the lone pairs of the 
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acceptor, 8 c i w 0, 9 a » /3 /2 = 54.74° in our notation. 
Actual angles deviate slightly from these ideal values and 
are known with some scatter. 

For the experimental geometry we get Ep = —7.046 
kcal/mol. Adding 1.820 kcal/mol for the exchange and 
dispersion energy [VRT(ASP-W)III 8 value] we get equi- 
librium binding energy D e = 5.110 kcal/mol. The model 
also yields the total dipole moment of the dimer in ex- 
cellent agreement with experiment (Table I). 

Since Ep is only a part of the total interaction, which 
also contains exchange and dispersion terms, we fix Roo 
and analyze the orientation dependence (Fig. 2). The 
minimum is achieved at 9d = 2.14° and 9 a = 66.26°, 
which is close to, but should not be confused with the 
equilibrium hydrogen-bonded configuration, since other 
terms may shift the minimum. Rotation of either the 
donor by Af^ w —(3, or the acceptor by A8 a » — j3 
produces an alternative hydrogen bonded arrangement 
sketched under the local minima in Fig. 2. 

Table I. Equilibrium binding energy D e (kcal/mol) and 
dipole moment p dlm (Debye) for water dimer. [i± m is the 
component of /i dlm normal to the principal axis. ' Geometry 
is fixed at experimental values; "^projection on the principal 
axis. 7 
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In order to further assess the quality of the model, we 
analyze the energy variation along a path where the ex- 
change and dispersion terms vary little. We choose to 
rotate the donor by an angle <f> around the bridging OH 
bond (Fig. 2). Only a single proton then changes its po- 
sition and stays far from all the nuclei of the acceptor at 
all cj). 

Figure 3 shows excellent agreement with all three best 
pair potentials. Note the small (< 1 kcal/mol) total 
amplitude of the variation, which is not described by a 
simple cos <f> function. The overall agreement in the full 
range of <fi is better with the ab-initio-b&sed SAPT-5s 9 
potential (the inset). However, at small <j> we get a near 
coincidence with the other two curves, VRT(ASP-W)III 
and SAPT-5st, 10 which are both spectroscopically-tuned. 
This is not surprising, assuming the spectroscopic tuning 
is more sensitive to the region near the equilibrium. 

Explicit distributed polarizabilities (4) suggest an esti- 
mate of the dispersion energy. Due to the fast r~ 6 decay, 
the dispersion is dominated by two terms, Eqq oc oloolo 
and EflQ oc anao- Small an in the second term 
is compensated by the proximity of the donor hydro- 
gen to the oxygen of the acceptor. Neglecting disper- 
sion nonadditivity and assuming an universal scaling of 
the dispersion coefficient ~ zuaolb for A and B 
species, we get Eq Q = za /R 00 — 0.99 kcal/mol and 
Eho = \ za nOLo I (Roo ~ 6 = 0-40 kcal/mol for linear 



hydrogen bond. The total E D = 1.39 kcal/mol can be 
compared to 1.56 kcal/mol from Fig. 3 of Ref. 33. For 
this crude estimate we used z = 344 kcal/mol value for 
Ar. The factor | accounts for the anisotropy of an- 
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Fig. 3 Energy variation for water dimer with rotation of 
the donor around the bridging OH bond. 

Since the minimal model is constructed based solely 
on monomer properties, we may speculate that it should 
describe larger clusters as well, where the non-pairwise 
additivity of energy is important. 9 ' 29 Such nonadditivity 
results from self-consistency of all the induced moments 
in the cluster, 24 and may be relevant for the cooperativity 
of hydrogen bonding in protein secondary structures. 35 

This work makes a step towards a chemical model for 
polarization intermolecular forces by combining minimal 
atomic multipoles with distributed polarizabilities, which 
together yield a transparent model for polarization forces 
in water. Its success raises a question of broader applica- 
bility especially to polarization and hydrogen bonding in 
peptides and proteins, and in water-protein interactions. 
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